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1.0  SUMMARY 


Approximations  to  the  boundary  integral  equation  (BIE)  formulation  for 
acoustic  wave  propagation  permit  simulation  of  acoustic  waves  in 
layered  earth  models  with  three  dimensional  layer  boundaries.  The 
complete  BIE  solution  is  approximated  by  a  series  expansion  analogous 
to  the  more  familiar  generalized  ray  expansion  widely  used  in 
seismological  modeling.  A  layer  to  layer  propagation  algorithm  is 
presented  which  is  efficient  enough  to  perform  three  dimensional  wave 
propagation  on  a  modern  minicomputer  equipped  with  an  array 
processor.  With  an  efficient  propagation  algorithm,  iterative  methods 
for  computing  the  layer  coupling  are  feasible.  The  "ray  expansion" 
approach  is  most  useful  for  approximating  solutions  on  wave  propagation 
problems  in  which  multiple  interaction  between  boundaries  can  be 
ignored. 

The  approximate  BIE  method  is  applied  to  an  acoustic  model  of  a 
mountain  in  which  a  flat  layered  velocity  structure  is  overlain  by  three 
dimensional  topography.  For  the  solution  that  includes  primary 
reflection  from  the  layered  velocity  structure  and  their  corresponding 
interaction  with  the  topography,  amplitude  variations  between  several 
profiles  can  be  interpreted  as  they  relate  to  the  topography  along  the 
profile.  The  location  of  the  source  and  receivers  on  the  near  surface, 
low-velocity  material  introduces  strong  acoustic  waveguide  effects 
equivalent  to  Love  waves  in  an  elastic  medium.  Modeling  these  guided 
waves  requires  including  waves  that  reflect  from  the  subsurface  velocity 
structure  and  interact  with  the  free  surface  several  times.  These 
guided  waves  dominate  the  solution  over  the  source-receiver  geometry 
of  interest.  Peak  amplitudes  vary  by  a  factor  of  2  for  stations  spaced 
at  1  km;  apparently  the  result  of  subtle  changes  in  the  interference  of 
waves  that  have  interacted  with  the  free  surface  in  different  ways. 
The  effect  of  three  dimensional  structure  on  peak  amplitudes  and 
waveforms  of  these  guided  waves  will  be  difficult  to  predict. 
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2.0  INTRODUCTION 

Efforts  to  understand  the  detailed  structure  of  the  earth  requires  the 
capability  to  numerically  model  seismic  wave  propagation  in  complex  two- 
and  three-dimensional  structures.  Computational  techniques  for 
modeling  wave  propagation  in  complex  structures  range  from  simple  ray 
tracing  to  finite  difference  and  finite  element  methods.  Ray  tracing 
methods  are  elegant,  conceptually  simple  to  grasp,  and  form  the 
intuitive  basis  for  understanding  many  types  of  wave  phenomena.  As 
an  interpretation  and  modeling  tool,  they  are  indispensible  for 
understanding  the  results  of  more  complex  calculations.  For  more 

complete  solutions,  finite  difference  methods  have  received  widespread 
attention  (e.g.  Alterman  and  Karal,  1968;  Boore,  1972;  Kosloff  and 
Baysal,  1982,  Vidale,  1986).  These  methods  have  large  computational 
requirements  and  it  is  not  currently  feasible  to  treat  three-dimensional 
problems.  Further  advances  in  computer  technology,  particularly  in  the 
development  of  large  parallel  machines  will  remedy  this  situation. 
However,  alternative  approximate  methods  that  provide  insight  into 
complex  wave  phenomena  will  always  be  necessary. 

One  method  for  posing  wave  propagation  problems,  dating  from  the 
early  nineteenth  century,  is  the  Boundary  Integral  Equation  (BIE) 

technique  in  which  a  wavefield  in  a  region  of  interest  is  represented  by 
the  values  of  the  field  and  its  gradient  along  a  bounding  surface.  While 
finite  difference  methods  require  sampling  the  wavefield  in  all  spatial 
dimensions,  the  BIE  method  requires  that  the  wavefield  be  sampled  only 
along  surfaces,  which,  at  first  glance,  appears  to  be  a  huge  saving  of 
effort  in  a  three-dimensional  problem.  However,  while  the  interaction 
between  samples  in  a  finite-difference  model  is  local,  the  interaction 

between  samples  in  a  BIE  representation  is  global:  each  point  on  a 

surface  interacts  with  every  other  point  on  the  same  surface  and  the 
adjacent  surfaces.  Thus,  for  complex  interaction  of  several  boundaries, 
(i.e.  "lots  of  multiples"  in  raytracing  terminology)  it  is  not  yet  clear 
that  the  BIE  method  represents  an  advantage.  However,  for  problems 
in  which  the  interaction  between  boundaries  is  not  severe  (i.e.  if 
primaries  and  low  order  multiples  are  of  interest),  then  BIE  methods 
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are  an  attractive  alternative,  particularly  when  propagation  over  many 
wavelengths  in  a  homogenous  media  is  involved. 

Several  treatments  of  the  BIE  formulation  have  appeared  in  the 

geophysical  literature  in  recent  years.  Cole  (1980)  presented  a 
time-domain  formulation  for  two-dimensional  acoustic  problems.  For 
many  problems  of  seismological  interest,  this  approach  will  be 
unsatisfactory  since  it  is  difficult  to  suppress  edge  reflections  and 
include  realistic  attenuation.  A  frequency-domain  treatment  of  the 

two-dimensional  elastic  problem  was  studied  by  Ferguson  (1982).  Both 
methods  failed  to  handle  the  interaction  integrals  efficiently  enough  to 
extend  the  results  to  three  dimensions  with  current  computer 
technology.  Schuster  (1984)  has  studied  several  approaches  to 
efficiently  solving  the  BIE's.  He  presents  an  iterative  solution  that  is 
directly  analogous  to  familiar  ray  expansions  used  in  plane-layered 
media  and,  thus,  has  a  strong  intuitive  basis.  His  formulation  lends 
itself  to  a  variety  of  useful  approximations  and  hybridizations;  it  is 
virtually  identical  to  the  approach  that  we  have  adopted  in  the  following 
report. 

Apsel  et  al.  (1985,  1983)  has  formulated  a  frequency  domain  approach 
that  is  stable,  involves  no  matrix  inversions,  and  that  can  include  the 
effects  of  realistic  material  properties  such  as  attenuation.  The 
formulation  is  designed  to  calculate  the  total  response,  including  all 
kinematic  and  dynamic  effects  using  a  specially  designed  pertubation 
treatment.  The  method  appears  to  be  more  cost-effective  than  finite 
difference  algorithms,  however  computational  demands  for  a  three 
dimensional  model  are  still  much  too  large  for  routine  treatment. 

An  attractive  feature  of  the  BIE  method  is  that  it  readily  leads  to 
intuitively  appealing  approximations  that  allow  some  hope  for 
understanding  complex  wavefields  in  terms  of  simpler  wave  propagation 
concepts  such  as  rays,  diffraction,  reflection  and  transmission.  An 
example  is  the  Kirchhoff  approximation  (Hilterman,  1970;  Berryhill, 
1979;  Scott  and  Helmberger,  1982;  Mellman  et  al . ,  1982)  in  which  the 
BIE  equation  is  broken  up  into  separate  calculations  for  the  propagation 
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between  layer  interfaces  and  for  the  interaction  along  individual  layer 
boundaries.  The  self  interaction  is  approximated  by  treating  each  point 
as  if  the  incident  wavefield  is  a  plane  wave  and  the  interface  is  a  plane 
oriented  along  the  tangent  to  the  surface  at  that  point.  Interaction 
between  two  points  on  the  same  boundary  is  not  considered  and  so  this 
approach  is  only  valid  if  the  topography  along  an  interface  is  not 
severe  enough  to  introduce  multiple  scattering;  i.e.  waves  reflected 
from  one  part  of  an  interface  to  another  part  of  the  same  interface  are 
negligible.  For  the  models  that  we  will  treat  in  this  report,  we  will  use 
this  assumption. 

In  the  following,  we  will  first  review  the  BIE  formulation  as  presented 
in  Apsel  et  al .  (1983,  1985)  and  provide  details  of  some  useful 
approximations .  Then  we  will  apply  the  methods  to  some  realistic 
calculations  to  model  strong  ground  motions  on  a  mountain  where 
topography  is  expected  to  have  significant  effects  on  the  wave 
propagation . 
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The  boundary  integral  equations  describing  wave  propagation  through 
arbitrary  three-dimensional  elastic  multilayered  media  are  derived  in  two 
steps.  First,  wave  propagation  within  a  single  irregular,  homogeneous 
layer  is  described  by  integral  representations  using  the  full  space 
Green's  functions  with  properties  of  that  layer.  Second,  the  wavefields 
in  each  irregular  layer  are  constrained  to  interact  at  the  layer 
boundaries  to  satisfy  all  of  the  boundary  and  continuity  conditions, 
which  leads  to  a  system  of  integral  equations  for  the  unknown  boundary 
values.  Once  this  system  of  equations  is  solved  for  the  boundary 
values,  the  wavefield  may  be  calculated  at  all  of  the  receiver  positions 
of  interest  using  the  integral  representations  of  the  first  step. 

The  model  geometry  for  the  wave  propagation  problem  solved  in  the  BIE 
formulation  is  depicted  in  Figure  1  by  N  irregular  layers  overlying  a 
semi-infinite  half-space.  The  layers  are  allowed  to  pinchout  but  not  to 
cross  in  this  formulation.  Each  layer  is  characterized  by  constant 
shear  and  compressional  wave  velocities  and  constant  densities. 
Material  attenuation  may  be  introduced  by  allowing  the  velocities  to  be 
complex.  Wave  propagation  within  a  given  layer  is  expressed  in  terms  of 
the  Green's  functions  for  a  full-space  with  the  properties  of  that  layer. 
The  formulation  is  not  restricted  to  constant  material  properties  within  a 
given  layer,  although  the  Green's  functions  for  that  case  are  quite 
simple.  In  the  following,  the  formulation  is  presented  for  the  acoustic 


1.  Cross-sectional  model  geometry  for  layered  half-space  used 
in  BIE  method  formed  by  N  irregular  layers  overlying  a 
uniform  half-space,  with  each  layer  characterized  by 
constant  material  properties.  The  source  and  receiver  can 
be  located  anywhere  in  the  medium. 
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The  potential,  4*/  then  satisfies  the  wave  equation: 


V2^  -  _1_  9jj>  =  0 


(1) 


2  2 
a  9t 


where  a  is  the  acoustic  velocity  and  p  is  the  density. 


If  we  define  the  Fourier  transform  with  respect  to  time  as 


$  (x,w) 


00 

■/ 


+iiut 

<)>  (x,t)  e  dt 


-00 


then  we  have  the  Helmholtz  equation 


V2  *  +  k2  4»  =  0 


(2) 


2  2  2 
with  k  =  ut/a  . 


The  first  step  in  the  formulation  is  to  write  expressions  for  the 
potential  field  within  a  single  layer  without  consideration  of  the 
boundary  interaction.  In  layer  £  (£=1 ,2, . . .  ,N+1 ),  the  potential  must 
satisfy  the  homogeneous  (£  $  s)  or  inhomogeneous  (£=s)  equations  of 
motion  (depending  on  whether  or  not  the  source  is  located  in  layer  £) 
for  a  full-space  with  properties  of  layer  £.  The  Kirchhoff  integral 
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representation  (e.g.  Stratton,  1948)  or,  equivalently,  the 

Representation  Theorem  of  Elastodynamics  (deHoop,  1958)  provide 
an  expression  for  the  potential  located  anywhere  within 

volume  containing  layer  £  in  terms  of  integrals  of  the  potential  and 
its  normal  field  over  the  bounding  surface  of  volume  multiplied  by 
the  corresponding  Green's  functions  for  a  full-space  with  properties  of 
that  layer.  The  potential  at  location  x£  can  then  be  written  in  the 
frequency  domain  using  the  integral  representation  for  a  volume 
bounded  by  layer  interfaces  S£  and  S^+1 . 


e  (x£)0(x 


g*  /(y£) 


H  (xt,y,)  »  (yt) 


■/ 


£+1 


..  £  ,  -*  » .  1 ,  , 
H  (VV,,,)*  (V„1> 


dS<vJ-*1) 


f 


+  6£s  G  < 


dV(z£) 


(3) 


in  which  the  frequency  arguments  have  been  omitted  for  brevity  and 


ym  =  an  integration  point  on  bounding  surface  Sm; 
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£  4  4 

G  (x^,ym)  =  the  full-space  Green's  function  for  the  potential 

#  at  location  y  on  surface^  S  due  to  point  force  in 
the  i-direction  at  location  x.  with  properties  of 
layer  £;  * 


♦*  (-*m)  =  ;  <Vm> 

^  =  'v  G‘  'vv.>  •  5<v*> 


f  (z^)  =  the  source  function  at  location  z£  anywhere  in 

layer  £  (assuming  the  source  is  a  Delta-function 
in  space,  then  the  volume  integral  reduces  to  the 
evaluation  of  the  integrand  at  point  z£); 


0,  if  £=s 
1 ,  if  £=s 


,s  =  source  layer  number; 


i,-*  . 

E 


1 ,  if  x£  inside  layer  £ 

if  x£  on  surface  bounding  layer  £ 
0,  if  x£  outside  layer  £ 


In  Eq.  (1),  the  layer  comprising  volume  V£  is  assumed  to  extend  to 
infinity  at  the  horizontal  extremes  to  eliminate  the  surface  integrals 
along  those  portions  of  the  surface  bounding  volume  and  the 

negative  sign  for  the  integral  over  surface  Sp  +  1  is  associated  with  using 
the  upward  normal  v  =  -v  in  the  definition  of  the  traction 

4  ■“ * 

components.  Once  the  boundary  values  for  4>(yrn)  and  <Kvm)  are 
determined  for  bounding  interfaces  and  S£  +  1 ,  Eq.  (1)  can  then  be 
used  to  obtain  the  displacement  field  at  any  point  x£  within  layer  £. 
Expressions  for  the  full-space  elastic  Green's  functions  with  constant 


SGI  -R-87-133 


10 


t, 

I 


material  properties  are  given  in  Appendix  A  of  Apsel  et  al . ,  1983  for 
two  and  three  dimensional  wave  propagation  in  elastic  as  well  as 
acoustic  media.  This  completes  the  propagation  step  of  the  BIE 
formulation  and  it  remains  to  impose  the  boundary  interaction  coupling. 

The  boundary  interaction  coupling  requires  simultaneous  satisfaction  of 
a  pressure-free  surface  (interface  1)  and  both  continuous  pressure  and 
normal  components  of  velocity  across  each  layer  interface  (2,3, . . .  ,N+1 ). 
The  coupled  boundary  integral  equations  arising  from  the  continuity 
conditions  across  each  layer  interface,  S£  (£=2,3, . . .  ,N+1 ),  are  obtained 
by  evaluating  Eq.  (4)  in  volumes  and  V£  (layers  £-1  and  £)  at  a 

discrete  number  (q£)  of  observations  along  common  surface  and 
imposing  the  boundary  conditions: 


0-1  -4  0  ^  0  •  1 

0  (x£)  =  4>  (x£)  and  _1_  0 

p£-1 

for  all  quadrature  points  y£  on  surface  S^.  To  derive  an  integral 

representation  for  0  when  0  is  on  the  boundary  requires  special 

£ 

treatment  of  the  singularities  present  in  the  Green's  functions  G  and 
£ 

H  .  Integration  of  the  singularities  is  straightforward  and  is  discussed 
in  work  by  Scott  (1985),  Cole  (1980)  and  Banaugh  (1962). 

The  resulting  set  of  equations,  for  £  =  1  to  N+1,  are 


n 


'•  <«*>  ■ 


and 


J 


£-1 


Vm)* 


'£-1 


£-1  -*  -*•  £-1  -> 

H  (x£'v£-1)  *  (Vje_1> 


dS(>£-1> 


)  - 


H  <xry£)  *  (yt) 


dS(yt)  * 


6*-1sF<*£> 


(4) 


Jf 


G£  (xry£)/(y£) 


n£  ^£'V£)o£(y£) 


ds(y£) 


(5) 


■J[ 


G"  <X£'V£+1)^(V£+1) 


£+1 


Hii  <VW  ♦  <y,.,> 


dS(vf<|) 


♦  6£sF(x£)  '  ( i ,  j  =  1 ,2 , 3 ) ,  (£=2,3 _ N  +  1  ) 


J 


where  F(xc)  =  |  G*(*£»z£K(2£)  dV  (z£). 
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Equations  (4)  and  (5)  represent  a  simultaneous  set  of  2q^  Fredholm 

integral  equations  of  the  second  kind  for  the  same  number  of  unknown 

£  £ 

boundary  values,  <$>  and  41  on  surface  S^,  which  are  coupled  to  the 
unknown  boundary  values  on  surface  ^  and  S^^  through  the 
integrals  over  surfaces  S^^  and  ,  respectively.  Note  that  when 

£=1  or  £=2  in  Eq.  (4),  then  the  term  involving  0\y^)  is  identically  zero 
because  of  the  tractionless  free-surface  conditions.  Also,  note  that 
when  £=N+1  in  Eq.  (5),  then  the  integral  over  surface  S^^  vanishes  by 
virtue  of  the  radiation  conditions  implicit  in  the  Green's  functions  for 
the  underlying  semi-infinite  space.  Boundary  integral  equations  (4) 
and  (5)  are  evaluated  at  a  discrete  set  of  q£  example  points  on  each 
boundary  S£  with  sample  points  from  one  boundary  becoming  Green's 
function  quadrature  points  for  the  next  boundary.  When  completely 
discretized,  Equations  (4)  and  (5)  represent  a  coupled  system  of 
singular  Fredholm  integral  equations  of  the  second  kind  for  the 
unknown  boundary  values  along  all  layer  boundaries.  Singularities 
occur  in  the  Green's  functions  when  quadrature  point  y  approaches 
observation  point  xm  in  the  second  integral  in  Eq.  (4)  and  the  first 
integral  in  Eq.  (5).  Singularities  can  also  occur  in  the  Green's 
functions  when  two  adjacent  layer  boundaries  intersect  and  would  be 
evidenced  in  the  first  integral  in  Eq.  (4)  or  the  second  integral  of  Eq. 
(5),  depending  on  where  the  intersections  occur. 

To  compute  a  numerical  approximation  to  (5)  we  discretize  the  integral 
equations  to  form  matrix  operations.  First  it  is  convenient  to  introduce 
some  additional  shorthand: 


Gm. 

') 

i 

= 

Gm(x.,  y  u  (yj) 

dS(y.) 

Hm 

>1 

= 

Hm(  x.,y  )  0  (y.  ) 

dS(  y . ) 
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acoustic  properties  am,  p^,  (m  =  i  or  j).  The  forcing  term  can  be 

written  as 


ik 


F.. 

•J 


I 


f  (z})  C  (x.,z.)  dZj 


where  F..  is  the  response  at  point  x  on  boundary  i  due  to  a  distribution 

„  -►  .  th 

of  sources  on  z  in  the  j  layer. 


We  define  submatrices,  D  (downgoing)  and  U  (upgoing): 


££-1 


l£-1 


.£-1 

’££-1 

0 


U 


££+1 


H 


££+1 


-G 


££+1 


the  boundary  interaction  matrix: 

,£-1 

I  r 

B, 


££ 


-H 


££ 

£ 

££ 


-G 


£-1 

££ 

.£ 

*££ 


1/2 

1/2 


and  the  forcing  function  matrix: 

’£-1  £ 

£  £ 

where  I  is  the  identity  matrix. 

The  unknown  boundary  values  are  written: 
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Then  the  integral  equation  (5)  can  be  written  as  a  system  of  equations: 


Note  that  the  free  surface  boundary  conditions  can  be  included  by 
explicitly  constraining  4>j 

and  that  the  radiation  conditions  are  explicitly  included  by  not  allowing 
any  upgoing  waves  from  the  N+1  region. 

Despite  a  simple  form,  Equation  (6)  represents  a  formidable  computational 

effort.  Each  submatrix  is  dense  and  for  problems  of  seismological 

interest,  say  with  grids  of  1000  x  1000  nodes,  when  propagation  over 

many  wavelengths  is  of  interest,  the  U,  D,  and  B  matrices  are  of 
6  6 

order  10  x  10  in  size.  Obviously  this  formulation  is  only  of  symbolic 
utility  and  some  approximations  are  necessary.  Furthermore,  decomposition 
into  simpler,  although  approximate,  methods  is  useful  from  the  standpoint 
of  understanding  the  complex  wave  interactions  in  terms  of  the  physical 
parameters  of  the  model. 

Several  ways  to  approximate  a  solution  to  Equation  (6)  may  be  suggested 
by  considering  simple  iterative  solutions  to  ordinary  systems  of  equations 
of  the  same  form.  Formally,  these  solutions  are  known  as  Born 
approximations  (Schuster,  1983)  but  here  we  use  a  heuristic  approach 
to  develop  the  same  ideas. 
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Consider  the  well  known  Gauss-Seidel  iteration  method  for  ordinary 
linear  equations  where  the  i  n  iteration  is  given  by: 


££- 1 


(i-1) 
£- 


U 


££+1 


4> 


(i-1) 

£+1 


(7) 


o 

Now,  on  the  first  iteration,  <t> 
row  and  moving  downward, 

At  that  layer 


is  set  equal  to  0.  Beginning  on  the  first 
remains  0  for  £<s  where  s  contains  a  source. 


and 


(1) 


B 


ss 


(D 

s+1 


B 


-1 

S+1 s+1 


s+1 


D 


s+1  s 


4> 


(D 


and  so  on  through  the  layers: 


4> 


(D 

N+1 


B 


N+1  N+1 


D  4>  ^  ^ 

UN+1  N  N 


At  the  end  of  the  first  iteration  down  through  the  layers,  we  calculate 
the  second  iteration  from  the  N+1  interface  up  through  the  layers 


4> 


(2) 


B'1  U  *(2) 

££+1  £+1 


16 


This  formulation  can  be  recognized  as  very  similar  to  a  ray  expansion 

with  the  U  and  D  operators  propagating  the  wavefield  from  one 

-1 

interface  to  the  next  and  the  determining  the  transmitted  and 

reflected  wave  field  at  each  interface.  It  is  identical  to  the  Generalized 

Born  Series  approach  described  by  Schuster  (1984).  If  the  layer 

interfaces  are  horizontal,  the  wavefield  can  be  Fourier  transformed  into 

a  horizontal  wavenumber  spectra  and  the  convolution  operators  U,  D, 
-1 

and  B  can  be  easily  computed  and  applied  for  each  wavenumber.  The 

matrix  equation  then  reduces  to  the  familiar  propagator  formulation  as 

described  in  textbooks  on  seismology  (e.g.  Aki  and  Richards,  1980). 

In  general,  however,  application  of  each  of  these  operators  represents  a 

convolution  in  either  the  spatial  or  wavenumber  domains.  Each 

2 

convolution  requires  on  the  order  of  q  operations  where  q  is  the 

3 

number  grid  nodes  on  each  surface.  For  a  model  of  grid  size  10  x 
3  12 

10  ,  this  represents  10  operations.  Clearly  further  approximations  to 
each  of  the  operators  are  necessary. 
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4.0  APPROXIMATIONS  TO  D£_1£,  U££  +  1 

To  propagate  a  wavefield  <J>£  from  one  layer  boundary  to  an  adjacent 
layer  boundary,  we  will  use  a  Fourier  transform  method.  Consider  the 
situation  in  Figure  2  in  which  we  wish  to  propagate  4>£  to  layer 
boundary  £+1 .  If  the  layer  boundaries  £  and  £+1  were  horizontal  we 
could  apply  a  two  dimensional  Fourier  transform  and  propagate  the 
wavefield  by  simple  multiplication  in  the  wavenumber  domains. 

The  use  of  FFT's  to  compute  the  wave  number  spectra  of  the  wavefield 

allows  us  to  calculate  this  operation  in  order  n  log  n  operations,  where 

n  is  the  number  of  points  defining  the  boundary,  a  large  savings  over 
? 

the  explicit  n  operations  in  the  (x,  y,  z)  domain. 

Now  consider  the  surface  z£  (x,y),  a  single-valued  function  of  x  and  y 
and  4>£(x,y),  the  sampled  wavefield.  We  sample  <t>£(x,y)  along  the 
intersection  of  z£(x,y)  and  a  horizontal  plane  =  constant: 


4>£(x,y'^k) 


Z£ 

\  1  dxdy 

l 

f°r  H  =  1 

^  / 

(8) 

=  0  for  z£  otherwise. 


where  is  the  distance  between  planes  on  which  a  sampling  of  is 
desired  and  v^  is  the  z  component  of  the  surface  normal.  if 
(x,y,z£)  is  sampled  over  K  horizontal  planes  such  that  <  z£  <  4K 
then 


D£+£  *£  (X'V'Z£) 


(9) 
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Figure  2.  Cross-section  description  of  propagation  algorithm.  Field 
values  o  and  V4>-v  for  4^.-]  <  +  1  (shaded  regions;  are  projected 

onto  a  plane  4^  =  constant.  The  project  field  values  are  the 
propagated,  using  Fourier  transform  techniques  to  each  plane  (j=1,J) 
on  surface  z^fx,^)  and  projected  onto  z^. 


where  now  represents  the  operator  that  propagates  from 

horizontal  level  to  layer  S. . 

After  representing  4>£  (x,y,z^)  in  terms  of  field  values  on  horizontal 

planes,  we  propagate  these  field  values  to  horizontal  planes 

=  =  t-  +  At  and  use  these  field  values  to  estimate  the  actual  values 

£+1  * 

of  ♦  ^  on  2£+i  using  a  simple  linear  interpolation  method  that  is  the 

inverse  of  (8)  above. 

Clearly,  the  number  of  operations  needed  to  propagate  the  wavefield  to 
or  from  a  surface  is  directly  related  to  the  relief  on  the  layer.  Since 
the  sampling  and  interpolating  is  analogous  to  that  used  in  some  simple 
finite  difference  schemes,  a  sample  interval  of  about  10  levels  per 
wavelength  will  be  necessary.  Therefore  if  the  total  relief  is 
comparable  to  many  wavelengths,  this  technique  becomes  less 
advantageous.  However,  it  should  be  remembered  that  in  our  frequency 
domain  formulation,  the  sampling  can  be  adapted  to  the  frequency  of 
interest,  with  only  coarse  sampling  necessary  for  the  low  frequencies. 
Nevertheless  for  most  problems  of  seismological  interest,  variation  in  the 
vertical  direction  is  much  less  than  the  horizontal  dimensions;  handling 
wave  propagation  in  homogeneous  layers  using  this  technique  may  allow 
us  to  treat  a  limited  set  of  realistic  models  without  requiring  a  super¬ 
computer. 
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5.0  APPROXIMATION  TO 


The  operator  is  very  large  and  completely  dense  (no  non-zero 

elements).  Calculating  the  inverse  directly  is  not  feasible  in  three 

dimensions;  even  calculating  the  elements  and  applying  the  operator 

directly  is  barely  manageable.  We  will  apply  the  propagation  scheme 
described  in  the  last  section  to  compute  B£^  <J>^,  which  may  provide  an 
iterative  means  of  calculating  B^  4>£. 

Consider  a  field  incident  on  a  boundary  £,  from  the  £-1  layer,  with 
O  o 

field  values  $  ,  V<t>  *v.  Further  suppose  that  we  have  an  estimate  of 

8,  0  ->  0  -+  ■+ 

the  boundary  values  0  (x)  =  0  (x)  and  (x)*v.  An  initial  estimate 
may  be  <j>e(x)  =  0°(x).  Then  the  discretized  boundary  interaction 

equations  can  be  written 


1/2 

j.£ 

*i 

=  -g*:1 

'j 

(a  Vo®v).  + 

H^."1  (b  0e).  ♦  <J>° 

ij  'j  yi 

1/2 

J- 

<J>i 

=  g£ 

'J 

(a  V<p®v). 

(10) 

e  ->  e  -* 
where  (b<J>  ).  =  b(y.)  <J>  (y.),  and  the  function  a(y)  is  defined 

J  J  J  £ 

similarly.  A  similar  set  of  equations  is  implied  for  V<J».*v: 


C-*  -*  /  p - 1  0  e  \  o 

1  2  Vo  *v  =  vV  -G  (aV0e-v)  ♦  H  .  ( b<p  ).  +  V<J»  -v. 

i  \  ij  J  'J  I  /  »  ' 


(11) 


1/2  V4>k’-v 


v.-V  G  (aV$  -v). 
'  \  'J  I 


o  0 

h  (b<t>  ). 
')  J 
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We  propose  to  iteratively  solve  Equation  (TO)  by  assuming  (1)  that  the 
♦£  depends  largely  on  a(y)  and  b(y)  near  y  =  x.,  and  (2)  that  a(y)  and 
b(y)  are  slowly  varying  functions.  Then  we  can  write 


MO  J-  - 

1/2  *i  "  _a(i)  Si 

+ 

£- 

b..N  h 
(0  i 

1/2  ♦?  *  a(i)  gf 

- 

bo)  hf 

where  (i)  indicates  no  summation  is  implied  over  i. 
In  Equation  12,  g^,  h^1  are  defined  by 


g0  = 

GP. 

'J 

(Vo1 

h0  = 

Hn. 

o>e 

1 

'J 

] 

For  Equation  11  we  have 


1/2  V<(.£  -v.  = 

ri  i 


-a.  v.*V  q.  1 

i  i  3 1 


b  v  *V  h£_1 
i  i  i 


+  V<{> 


1/2  V0£  v. 

i  i 


a.  v.*Vg. 

i  i  3i 


b.  v.'V  h£ 

i  i  i 


The  key  assumption  in  the  above  manipulations  is  the  removal  of  the 
unknown  functions  a(y)  and  b(y)  from  the  convolution  of  the  Green's 


functions,  G..  and 
U 


H. 


•j 


with  the  estimated  field  values 


and  V<j>  *v. 


These  convolutions  can  now  be  computed  in  a  way  exactly  analogous 
to  the  methods  described  in  the  previous  section. 
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Now  the  continuity  conditions  at  the  boundary  give  two  equations  for 
the  unknowns  a.  and  b.  at  each  grid  point,  which  will  define  the 
functions  a(x)  and  b(x)  necessary  to  update  estimates  of  the  field 
values . 

S.  -* 

By  applying  the  boundary  conditions  to  <f>.  and  V<j>.  *v.,  we  can  solve  for 
the  function  a(x),  b(x)  for  all  surface  points  providing  a  new  estimate 

<jie  (x)  =  a(x)  4>e  (x)  and  V<t>e  (x)  •  v  =  V<pe  (x)  •  v 

new  old  new  old 

The  procedure  is  then  iterated  several  times.  The  propagation  methods 
discussed  in  the  previous  section  to  propagate  field  values  from  layer  to 
layer  can  be  extended  to  propagate  field  values  between  different  areas 
of  the  same  interface.  The  above  iterative  solution  has  undergone  only 
preliminary  exploration  and  thus  its  convergence  properties  are 
currently  unknown.  It  is  an  appropriate  way  to  start  investigating 

methods  to  solve  the  boundary  interaction  problem  given  a 

computationally  feasible  way  to  apply  the  operator  B££. 

The  models  provided  by  AFGL  that  we  will  use  in  the  subsequent 

-1 

sections  of  this  report  require  no  elaborate  approximations  to  B^.  The 
only  layer  boundary  that  varies  in  these  dimensions  is  the  free  surface 
boundary.  Furthermore,  the  topography  is  smooth  enough  to  permit 
some  reasonable  approximations  that  represent  huge  saving  in  the 
calculations.  At  the  free  surface,  the  boundary  condition  is  4>=o.  We 
make  the  further  assumption  that  V0^’V  =  -V0O,v,  where  the 

superscript,  R,  designates  the  reflected  field  and  o,  the  incident  field. 
This  assumption  neglects  energy  that  may  reflect  several  times  off  of 
the  surface  during  a  single  interaction,  such  as  in  the  "whispering 

gallery"  effect.  This  approximation  is  identical  to  the 

Kirchhoff-Helmholtz  approximation  as  described  by  Scott  and  Helmberger 
(1983). 
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6.0  AFGL  MODEL  CALCULATIONS  j 

The  earth  structures  used  for  the  wave  propagation  studies  are  from  | 

Cipar  (Personal  Comm.)  and  referred  to  as  Generic  Mountain,  model  1 
and  model  2.  Table  1  gives  the  model  parameters,  the  S-wave  velocity 
being  irrelevant  for  our  purposes  because  we  are  dealing  with  acoustic 
waves  (i.e.,  we  assume  that  the  model  is  a  fluid).  The  three 
dimensional  variation  of  the  model  is  confined  to  the  free  surface;  the 
lower  layers  are  flat.  As  discussed  earlier,  this  feature  greatly 
simplifies  the  calculations  since  propagation  and  interaction  between  flat 
layers  can  be  handled  with  Fourier  transforms. 

Figure  3  gives  a  contour  map  of  the  surface,  and  Figures  4,  5,  and  6 
show  several  cross-sections  through  the  model.  The  actual  model  used 
in  the  calculation  is  a  smoothed  version  of  that  provided  b^  Cipar 
(Personal  Comm.).  The  contours  were  digitized,  and  converted  to  a 
regularly  sampled  grid  using  a  surface  inversion  algorithm  that 
determines  the  smoothest  unaliased  surface  that  fits  the  data  values. 

This  gridding  method  removes  the  sharp  corners  found  in  the  model  as 
described  by  Cipar  (Personal  Comm.)  and  creates  a  model  more  suitable 
for  numerical  calculations. 

The  source  locations  and  receiver  arrays  are  shown  in  Figure  4  and 
described  in  Table  2.  In  addition  to  the  receiver  arrays  specified  by 
AFGL  we  have  computed  several  profiles  at  the  receiver  depth  in  the 
directions  of  the  receivers  arrays.  These  profiles  are  useful  for 
identifying  the  origin  of  the  arrivals  based  upon  the  moveout  of  the 
various  phases. 

The  maximum  feasible  grid  size  for  the  calculation  is  currently  128  \ 

128.  To  obtain  the  maximum  frequency  bandwidth  w<e  performed  the 
calculations  using  two  different  grid  node  spacing  depending  upon  the 
frequency.  A  single  grid  with  node  spacing  of  0.24  km  was  used  to 
calculate  the  response  at  the  three  profiles  from  0  to  5.4  hz.  For 
higher  frequencies,  5.6  to  11  hz,  separate  grids  were  used;  one  tor  the 
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TABLE  1 


MODEL  1 


SEISMIC  VELOCITY 

DENSITY 

LAYER 

P-WAVE 

S-WAVE 

g/cm3 

1 

2.6 

1.49 

2.65 

mountain 

2 

3.2 

1.86 

2.65 

basement 

3 

6.0 

3.  SO 

2.80 

granite 

MODEL  2 

1 

2.60 

1.49 

2.65 

cap rock 

2 

2.21 

1.27 

2.25 

low  velocity 

3 

3.20 

1.86 

2.65 

basement 

4 

6.00 

3.50 

2.80 

granite 
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EAST-WEST  DISTANCE  IN  KILOMETERS 


FIGURE  3.  Contour  map  of  surface  topography  of  generic  mountain 
model.  Depth  units  are  in  km. 
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TABLE  2 


SOURCE  -  RECEIVER  GEOMETRY 


X  (N-S)  Y  (E-W)  Z  (DEPTH)  LAYER 

(km)  (km)  (km) 

Source  25.0  20.0  0.8  1 


RECEIVER  PROFILES  AT  Z  (DEPTH)  =  0.4572  km 


Profile 

AZIMUTH* 

DISTANCE  TO 

1st  RECEIVER 

DISTANCE  TO 
LAST  RECEIVER 

RECEIVER 

SPACING 

(degrees) 

(km) 

(km) 

( km) 

A 

338 

5 

10 

1 

B 

67 

3 

10 

1 

C 

116 

5 

10 

1 

ARRAY 

RECEIVER 

ARRAYS  AT 

Z  (DEPTH)  = 

0.4572  km 

AZIMUTH 

STATION 

X  (N-S) 

Y  (N-S) 

DISTANCE 

FROM  SOURCE 

(km) 

(km) 

(km) 

( degrees ) 

A 

1 

16 

16 

9.85 

336 

A 

2 

16 

17 

9.49 

342 

A 

3 

17 

16 

8.94 

333 

A 

4 

17 

17 

8.54 

339 

B 

1 

22 

25 

5.83 

59 

B 

2 

22 

26 

6.71 

63 

B 

3 

23 

25 

5.39 

68 

B 

4 

23 

26 

6.32 

72 

C 

1 

28 

27 

7.62 

113 

C 

2 

28 

28 

8.54 

in 

C 

3 

29 

27 

8.06 

119 

C 

4 

29 

28 

8.94 

117 

♦Source 

receiver  azimuth  measured 

clockwise  from 

N . 
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A  profile  and  one  for  the  B  and  C  profiles,  with  1/2  the  grid  node 
spacing  used  in  the  low  frequency  calculation.  The  high  and  low 
frequencies  for  each  receiver  response  were  then  merged.  For  each 
calculation,  the  model  was  sampled  on  a  128  x  128  grid,  centered  on  the 
area  of  interest.  The  grid  was  smoothly  tapered  to  a  horizontal 
boundary  on  the  edges  and  extended  to  256  x  256  points  for  the 
calculation.  The  frequency  range  of  the  calculations  were  from  0  to  11 
hz  at  an  interval  of  0.2  hz  giving  a  5  s  time  window  for  each  response. 

Tapering  and  padding  the  model  suppressed  only  some  of  the 
wraparound  and  edge  effects  resulting  from  using  the  Finite  Fourier 
Transform  in  the  calculations.  We  are  still  faced  with  aliasing  and 

wraparound  effects  that  cannot  be  removed  by  simply  increasing  the 
size  of  the  model.  These  result  from  the  implicit  spatial  and  temporal 
periodicity  of  the  wavefield  in  the  frequency  domain  using  finite  Fourier 
transforms:  the  source  is  assumed  to  be  both  periodic  in  time  and 

located  at  regular  intervals  in  space.  The  amplitude  decay  with 
distance  between  source  and  receiver  that  results  from  the  wave 

equation  cannot  sufficiently  attenuate  arrivals  from  adjacent  intervals 
and  from  sources  other  than  the  one  of  interest.  The  numerical 

techniques  necessary  to  describe  three  dimensional  wave  propagation 
using  Fast  Fourier  Transforms  is  described  by  Bouchon  (1979)  and  is 
further  discussed  in  Aki  and  Richards  (1980).  Essentially,  the 
response  is  computed  using  a  complex  value  for  frequency,  of  which  the 
imaginary  part  is  assigned  to  suppress  the  amplitude  of  propagating 

waves  with  distance.  After  computing  the  time  domain  response,  the 
effects  of  the  complex  frequency  value  can  be  removed  for  the  time 
window  of  interest  by  multiplying  by  an  exponentially  increasing  time 
function . 
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7.0  RAVTRACING  AND  BIE  CALCULATIONS  FOR  AFGL  MODEL 


To  check  the  results,  we  compared  the  profiles  with  results  from 

TM 

Sierra's  raytracing  package  QUIKSHOT  The  BIE  formulation 

presented  earlier  makes  these  comparisons  easy  since  we  can  effectively 
compute  a  ray  expansion  and  explicitly  include  the  same  propagation 
paths  modeled  in  the  raytracing.  To  trace  rays  we  used  model  2  as 
described  in  the  previous  section  and,  again,  simulated  a  fluid  medium 
by  making  the  S-wave  velocity  negligibly  small.  We  designed  the 
raytracing  problem  to  be  simple  enough  to  easily  interpret  the  results, 
yet  complex  enough  so  that  complications  due  to  three  dimensional 
effects  can  be  checked  and  compared.  We  compute  what  we  will  call  the 
primary  response,  which  includes  all  primary  reflections  from  the 
subsurface  layers  as  well  as  the  direct  wave.  Because  of  the  proximity 
of  the  source  and  receivers  to  the  free  surface  we  include  all  of  the 
interaction  with  the  free  surface  near  both  the  source  and  the 
receivers.  In  all,  13  raypaths  were  deemed  necessary  to  compute  a 
"primary"  solution;  these  are  schematically  given  in  Figure  7. 

The  comparison  of  the  BIE  and  raytracing  results  are  presented  in 
Figures  8  and  9.  The  source  function  in  all  cases  is  a  "Ricker" 
wavelet,  shown  in  Figure  10,  which  is  the  second  time  derivative  of  a 
Gaussian  pulse  with  a  peak  frequency  response  at  5.5  hz.  We  have 
plotted  both  pressure  and  the  vertical  component  of  displacement. 
There  is  some  disagreement  between  the  results  at  receivers  6,  7  and 
8;  note  the  larger  relative  amplitude  of  the  second  arrival.  This  is 
because  the  rays  interacting  with  the  lower  layers  go  through  the 
critical  angle  at  these  distances.  The  amplitudes  of  the  raytracing 
results  are  least  accurate  near  the  critical  angle.  Overall,  the 
agreement  between  the  raytracing  solution  and  the  BIE  solution  is  very 
good,  providing  confidence  that  the  BIE  calculations  are  working 
correctly.  Conversely,  it  shows  that  ray  tracing  provides  very  useful 
results  for  some  problems. 

To  compare  the  AFGL  models,  Model  1  and  2  from  Generic  Mountain  B, 
we  show  the  three  profiles  for  the  two  models  in  Figures  11  to  16. 
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BIE  RESULTS:  PROFILE  B.  Z  COMPONENT 

HORIZONTAL  SEISMIC  SECTION  AT  DEPTH  (  1)=  0.45720 

TYPE  2:  X.  Y  (151) «  23.  91  22.  73  T0X.Y(158)=  21.09  i 


FIGURE  8. b. 


RflYTRflCING  RESULTS:  PROFILE  B 

4  LAYER  GENERIC  MOUNTAIN  MOOEL 
RESSURE  COMPONENT.  5.5  HZ  RICKER  WAVELET 


FIGURF  9. a.  Comparison  of  pressure  component  synthetic  seismograms 
computed  by  (a)  raytracing  and  (b)  the  Bit  method  for  profile  B. 


BIE  RESULTS:  PROFILE  B.  PRESSURE  COMPONENT 

HORIZONTAL  SEISHIC  SECTION  AT  DEPTH  I  I)  -  0.45720 
TYPE  2:  X .  Y  (  25)-  23.  91  22.73  TO  X.Y(  32)-  21.09  29.06 


FIGURE  9.b. 
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BIE  RESULTS:  MODEL  1.  PROFILE  R.  RADI RL  COMPONENT 

HORIZONTAL  SEISMIC  SECTION  NT  OEPTH  (  II-  0.45720 
TTPE  2s  X.TI  59)  •  20.39  18.05  TO  X.TI  641-  IS.  70  16.17 
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BIE  RESULTS:  MOOEL  2.  PROFILE  R.  RRDIRL  COMPONENT 

HORIZONTAL  SEISMIC  SECTION  NT  OEPTHI  II-  0.45720 
TTPE  2«  X.TI  591-  20.39  19.05  T0X.TI84I-  15.70  19.17 


FIGURE  12.  The  radial  component  of  the  displacement  response  for 

model  2,  profile  A,  calculated  by  the  BIE  method  for  primary  reflections 

-14 

only.  The  traces  are  scaled  to  a  peak  amplitude  of  3.63  x  10  m  for 
trace  1 . 
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These  calculations  were  made  for  the  primary  response,  which  considers 
only  those  raypaths  (or  energy  paths)  diagrammed  in  Figure  7.  In  all 
cases  we  give  the  radial  component  of  displacement,  since,  for  the 
geometry  under  consideration ,  this  will  be  the  dominant  component  of 
motion.  The  amplitudes  are  scaled  to  a  pressure  source  with  a  step 
function  time  history  of  1  N-m  amplitude.  The  seismograms  have  been 
low-pass  filtered  with  a  frequency  taper  between  8.8  and  11  hz.  Each 
profile  extends  from  about  5  to  10  km;  the  plane  of  the  receivers  (at  Z 
=  0.4572  km)  intersected  the  free  surface  near  4  km  along  profiles  A 
and  C.  Therefore,  the  closest  receiver  located  below  the  surface  of  the 
mountain  was  beyond  the  distance  where  critical  reflection  from  the  3.2 
km  interface  occurs  for  both  models  (this  critical  distance  is  2.5  km  for 
model  1  and  2.2  km  for  model  2). 

The  difference  between  the  two  models  (i.e.,  the  low  velocity  zone) 
shows  up  in  differences  in  the  peak  amplitude  of  the  first  arrivals.  In 
model  1  the  first  arrival  is  an  interference  o'f  the  direct  ray,  head  wave 
from  the  3.2  km  interface,  the  post-critical  reflection  from  the  3.2  km 
interface  and  all  of  the  associated  interactions  with  the  free  surface. 
These  arrivals  are  all  within  0.12  s.  The  firs  large  arrival  on  all  of 
the  model  1  profiles  is  the  result  of  the  interference  of  the  direct  wave 
and  the  post  critical  reflection  from  the  3.2  km  interface.  A  small  head 
wave  from  the  3.2  km/s  layer  becomes  evident  at  farther  ranges.  In 
model  2,  the  direct  wave  arrives  almost  0.2  s  before  the  post  critical 
reflection  from  the  3.2  km/s  layer,  which  results  in  a  more  dispersed 
group  of  first  arrivals  with  lower  peak  amplitudes  than  for  model  1.  In 
both  models  the  amplitude  of  the  later  reflection  from  the  6  km/s 
interface  is  similar;  this  arrival  has  a  critical  distance  of  about  8  km. 

While  the  difference  in  peak  amplitudes  between  models  1  and  2  can  be 
explained  by  the  presence  of  the  low  velocity  zone,  the  difference  in 
amplitudes  among  the  three  profiles  may  be  explained  by  the  structure 
of  the  free  surface.  One  large  difference  is  at  the  near  distance  range 
of  5  and  6  km:  profile  C  shows  markedly  smaller  peak  amplitudes  than 
profiles  A  and  B,  which  are  similar.  Reference  to  Figures  3,  4,  5  and 
6  show  that  two  effects  are  likely  important.  First,  for  profiles  A  and 


C,  the  receiver  at  5  km  distance  is  in  a  "shadow"  for  the  direct  wave, 
i.e.  a  ray  from  the  source  to  the  5  km  receiver  intersects  the  free 
surface,  which  reduces  the  amplitude  of  the  direct  arrival.  This  effect 
likely  explains  the  slight  amplitude  differences  between  profiles  A  and 
B.  Second,  profiles  A  and  B  encounter  largely  two-dimensional 
structure  (see  Figure  3)  for  the  wavelengths  of  interest;  in  other 
words,  the  energy  propagates  along  the  same  direction  as  the  major 
variations  in  the  structure.  The  free  surface  interaction  along  profile 
C  is  more  three  dimensional,  the  structure  varying  obliquely  to  the 
direction  of  wave  propagation,  and  energy  is  scattered  out  of  the  plane 
of  propagation,  significantly  reducing  the  amplitudes  at  the  near 
ranges. 

The  amplitude  pattern  among  the  profiles  changes  at  the  more  distant 
ranges  (9-10  km)  where  profile  8  shows  peak  amplitudes  less  than  1/2 
those  of  profile  A  and  C.  The  reason  is  similar  to  that  discussed 
above.  Profiles  A  and  C  are  encountering  two-dimensional  structure  at 
these  ranges  while  profile  B  encounters  three  dimensional  structure  that 
causes  significant  out-of-plane  scattering  in  the  free  surface 
interaction. 

Consideration  of  the  more  familiar  wave  propagation  phenomena  in  plane 
layered  media  indicates  that  at  the  distance  ranges  of  interest  in  this 
study,  the  primary  reflections  are  not  the  dominant  contribution  to  the 
seismogram.  The  distances  are  well  beyond  critical,  and  the  source  is 
located  in  a  low  velocity  layer,  which  means  that  considerable  energy  is 
propagating  as  trapped  waves  in  the  near  surface  layers.  In  our 
acoustic  model,  these  trapped  waves  are  analogous  to  the  Love  waves  in 
elastic  models  for  the  earth.  To  model  these  waves  the  calculations 
must  account  for  several  interactions  between  the  lower  layers  and  the 
free  surface,  requiring  considerably  more  computer  effort.  The 
convergence  of  the  ray  expansion  has  been  considered  by  several 
authors  (e.g.  Hron,  1971);  we  simply  argue  that,  since  10  km  is 
approximately  3  times  the  critical  distance,  we  can  reasonably  expect  to 
model  these  waves  by  considering  up  to  three  reflections  from  the  loner 
layers. 
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We  performed  the  more  complete  calculations  only  for  Model  2  because  of 
the  effort  involved.  Large  differences  between  the  two  results  are  not 
expected,  and,  given  the  level  of  approximation  involved  by  neglecting 
many  important  elastic  effects  (such  as  Rayleigh  waves),  it  is  not 
expected  that  meaningful  comparisons  between  the  two  models  would 
result.  Figures  17  through  19  show  the  radial  displacement  components 
for  the  three  profiles  for  Model  2.  In  all  instances  the  trapped  waves 
dominate  the  responses;  there  is  very  little  amplitude  decay  with 
distance.  Now  it  is  profile  B  that  shows  consistently  larger  amplitudes 
out  to  a  distance  of  9  km.  Again,  the  waveforms  for  profiles  A  and  C 
are  remarkably  alike,  although  the  amplitudes  on  the  near  traces  are 
about  25%  larger  on  profile  A.  The  three  dimensional  structure  between 
the  source  and  the  near  receivers  on  profile  C  is  likely  causing  the 
reduced  amplitudes;  the  similarity  in  waveforms  between  A  and  C  is 
remarkable  in  view  of  the  very  different  waveforms  for  profile  B. 

From  the  prefiles  presented  in  Figures  17-19  we  see  that  interpreting 
trapped  waves,  or  equivalently,  surface  waves  that  have  propagated 
through  arbitrary  three-dimensional  structures  will  be  difficult.  More 
modeling  studies  will  be  necessary  to  determine  data  analysis  techniques 
that  can  map  observable  features  onto  structural  parameters.  Faster 
computers  and  improved  algorithms  for  quickly  calculating  partial 
solutions  are  essential  for  learning  how  to  interpret  and  predict 
three-dimensional  wave  propagation  effects. 

To  verify  the  completeness  of  the  solution  within  the  context  of  an 
acoustic  model,  we  computed  complete  synthetics  using  a  wavenumber 
integration  method  (VESPA™)  developed  by  Apsel  (1979)  for  a 
flat-layered  acoustic  approximation  to  Model  2.  The  model  was  modified 
so  that  the  free  surface  was  at  a  depth  of  0.4  km,  with  the  source, 
receivers,  and  deeper  structure  unchanged.  The  results  are  shown  in 
Figure  20.  Although  a  seismogram  by  seismogram  comparison  is  not  very 


revealing,  many  of  the  features  in  the  profiles  computed  by  the  BIE 
method  are  found  in  the  VESPA  solution.  In  particular,  the 
approximate,  relative  decay  of  amplitudes  is  similar,  as  well  as  the 
duration  of  the  high  amplitude  portions  of  the  seismograms. 


SGI  -  R  -87- 1 33 


RESULTS:  MODEL  2.  PROFILE  B,  RflDIRL  DISPLACEMENT 

HORIZONTAL  SEISMIC  SECTION  AT  DEPTH  (  1)=  0.45720 

!:  X.  Y  (  67)  =  23.  91  22.73  TO  X.Y(  74)=  21.09  29.  06 
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18.  The  radial  component  of  the  displacement  response  for 
profile  B  calculated  by  the  BIE  method  including  multiples  to 
irder.  The  traces  are  scaled  to  a  peak  amplitude  of  5.29  x 
for  trace  2. 
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E  20.  The  radial  component  of  displacement  for  a  flat  layered 
calculated  using  a  wavenumber  integration  method  (Apsel,  1979)/ 
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We  present  the  calculations  for  each  of  the  three  receiver  arrays  for 
both  the  primary  and  extended  response  calculated  for  Model  2.  The 
traces  have  been  convolved  with  a  Ricker  wavelet,  peaked  at  5.5  hz,  to 
simulate  the  source  spectrum  as  suggested  by  AFGL;  the  wavelet  and 
its  spectrum  are  shown  in  Figure  10.  Figure  21  through  23  show  the 
vertical  and  radial  primary  responses  for  the  three  arrays.  Table  2 
gives  the  array  locations.  The  variations,  within  the  arrays,  of  the 
relative  amplitudes  and  the  waveforms  are  slight  indicating  that  the 
primary  response  is  stable  over  1  km  distances.  There  is  some 
amplification  at  station  B-3  over  what  would  be  expected  from  normal 
decay  with  distance  from  the  source.  Array  B  shows  the  largest  peak 
amplitude  variation,  only  part  of  which  is  due  to  source-receiver 
distance. 

The  amplitude  and  waveform  stability  of  the  primary  response  is  in 
contrast  to  the  variations  observed  in  the  extended  response.  Figures 
24-26  give  the  radial  and  vertical  displacement  response,  including 
multiples  up  to  the  third  order,  for  each  array.  Now  peak  amplitude 
variations  of  a  factor  of  2  between  stations  become  evident  for  arrays  A 
and  C.  This  variation  is  larger  than  that  associated  with  normal  decay 
of  amplitude  with  distance  as  can  be  seen  by  examining  Figure  20.  The 
amplitude  variations  are  the  result  of  subtle  differences  in  the  way  the 
multiple  responses  combine.  There  is  no  obvious  way  to  predict  these 
amplitude  variations  short  of  modeling;  it  is  very  difficult  to  explain  the 
behavior  in  the  modeling  results  in  terms  of  simple  wave  phenomena. 

For  the  frequencies  and  wavelengths  considered  here,  (<10  Hz,  >0.25 
km)  the  effects  of  surface  topography  on  the  amplitude  and  duration  of 
strong  shaking  is  not  large  compared  to  typical  uncertainties  in 
seismology  arising  from  unknown  velocity  structure.  Given  these 
unknowns,  the  amplitude  prediction  from  a  flat  layered  wave-number 
integration  method  or  even  a  generalized  ray  method  are  adequate. 
Given  the  dominance  of  the  waveguide  effects,  either  of  these  methods 
will  be  much  more  satisfactory  than  simple  ray  tracing. 


BIE  RESULTS:  MODEL  2.  ARRAY  A,  5.5  HZ  WAVELET 

HORIZONTAL  SEISMIC  SECTION  AT  DEPTH!  1)  =  0.45720 
TTPE  3:  X.T!  43)=  15.  94  15.94  TO  X.  Y  C130I  =  17.11  17.11 
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RESULTS:  MODEL  2,  ARRAY  C,  5.5  HZ  WAVELET 

HORIZONTAL  SEISMIC  SECTION  AT  DEPTH!  1)=  0.45720 
i  X.T(  51)=  27.89  26.95  TO  X.Y!138)  =  29.  08  27.  89 


GROUP  NUMBER 

FIGURE  23.  The  radial  and  vertical  components  of  displacement 
response  for  model  2,  profile  C  calculated  by  the  BIE  method  using 
primaries  only.  See  caption,  Figure  21.  The  traces  are  scaled  to  a 
maximum  amplitude  of  9.96  x  10  m  on  trace  1  (station  C-1,  radial). 


BIE  RESULTS:  M0OEL  2,  ARRAY  A,  5.5  HZ  WAVELET 

HORIZONTAL  SEISHJC  SECTION  RT  DEPTH  (  t)«  0.45720 

TYPE  3*  X.Yt  431-  15.94  15.94  T0X.Y1130)-  17.11  17.11 
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BIE  RESULTS:  MODEL  2,  ARRRT  B,  5.5  HZ  HflVELET 

HORIZONTAL  SEISMIC  SECTION  RT  DEPTH  I  1)=  0.  45720 
TYPE  3:  X.  Y  (  47)=  22.  03  25.08  TO  X.YU34U  22.  97  26.  02 


GROUP  HUMBER 

FIGURE  25.  The  radial  and  vertical  displacement  response  for  model 
2,  array  B  calculated  by  the  BIE  method  including  multiples  to  third 
order.  See  caption,  Figure  24.  The  traces  are  scaled  to  a  maximum 
amplitude  of  2.50  x  10  ^m  on  trace  1  (station  B-1,  radial). 
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8.0  CONCLUSION  AND  RECOMMENDATIONS 

The  modeling  results  presented  in  this  report  represent  an  initial  effort 
at  calculating  the  response  of  three  dimensional  structures.  We  have 
successfully  demonstrated  that  wave  propagation  in  a  limited  class  of 
acoustic  models  can  be  numerically  simulated.  Because  the  results  in 
this  report  rely  on  acoustic  models,  the  amplitude  predictions  for  what 
are  essentially  surface  waves  must  be  interpreted  carefully.  It  is 
expected  that  the  Rayleigh  waves  will  contribute  to  the  solution  as  much 
or  more  than  the  trapped  waves  in  the  above  calculations.  Since  these 
propagation  modes  travel  at  less  than  the  S-wave  velocity,  the  duration 
of  large  amplitude  motion  will  be  increased.  The  Rayleigh  waves  will 
also  be  most  affected  by  the  surface  topography  and  nothing  in  these 
calculations  is  useful  to  predict  these  effects. 

The  techniques  developed  in  this  report  will  be  useful  for  certain 
elastic  problems  in  which  multiple  interaction  between  boundaries  are 
not  important.  Problems  concerning  three  dimensional  modeling  of 
reflection  seismograms  of  teleseismic  body  waves  are  examples.  The 
propagation  algorithm  presented  in  this  report  has  proven  to  be 
economical  for  propagating  a  wavefield  from  one  boundary  to  another  in 
a  three  dimensional  structure.  It  remains  to  develop  an  algorithm  to 
efficiently  compute  the  boundary  interaction  terms. 
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